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Abstract 

A strong-coupling expansion for the Green's functions, self-energies and correlation functions 
of the Bose Hubbard model is developed. We illustrate the general formalism, which includes 
all possible inhomogeneous effects in the formalism, such as disorder, or a trap potential, as well 
as effects of thermal excitations. The expansion is then employed to calculate the momentum 
distribution of the bosons in the Mott phase for an infinite homogeneous periodic system at zero 
temperature through third-order in the hopping. By using scaling theory for the critical behavior 
at zero momentum and at the critical value of the hopping for the Mott insulator to superfluid 
transition along with a generalization of the RPA-like form for the momentum distribution, we 
are able to extrapolate the series to infinite order and produce very accurate quantitative results 
for the momentum distribution in a simple functional form for one, two, and three dimensions; 
the accuracy is better in higher dimensions and is on the order of a few percent relative error 
everywhere except close to the critical value of the hopping divided by the on-site repulsion. In 
addition, we find simple phenomenological expressions for the Mott phase lobes in two and three 
dimensions which are much more accurate than the truncated strong-coupling expansions and any 
other analytic approximation we are aware of. The strong-coupling expansions and scaling theory 
results are benchmarked against numerically exact QMC simulations in two and three dimensions 
and against DMRG calculations in one dimension. These analytic expressions will be useful for 
quick comparison of experimental results to theory and in many cases can bypass the need for 
expensive numerical simulations. 

PACS numbers: 03.75.Lm, 37.10.Jk, 67.85.Hj 
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I. INTRODUCTION 



The Bose Hubbard model 



was heavily studied as a simple model for disordered 
superconductors [2]; subsequently it was demonstrated 4| that ultra-cold atoms trapped 
in optical lattices provide an alternate, and more controllable, experimental realization of 
it, sparking even more interest. One of the most useful tools for analyzing the states of 
cold atom systems is a time-of-fiight measurement of their momentum distribution when 
the lattice and trapping potentials are rapidly shut off and the atomic cloud is allowed to 
expand and then is imaged with absorption spectroscopy The time-of-fiight image, in 
the long-expansion-time limit, is directly related to the momentum distribution function of 
the atoms in the optical lattice before expansion 

aa. 

Even before cold atom systems were employed to examine Bose Hubbard model physics, 
the phase diagram of the model was accurately determined in a strong-coupling approxima- 
tion 8|, |9| (for a recent review of this early history, see Ref. llOl). This approach, which relied 
on expanding the properties in a perturbative series in the hopping, captured much of the 
behavior of the model, and when extrapolated via a scaling theory ansatz for the critical 
behavior at the tips of the Mott lobes 2|], proved to be as accurate as the quantum Monte 



Since then, the 



14,^ 



Carlo (QMC) simulations that had been performed at that time 
strong coupling perturbation theory has been pushed to higher order |13l . 
QMC simulations have improved dramatically in two 16| and three dimensions [T 
addition, highly accurate density matrix renormalization group (DMRG) studies have been 



performed on the model in one dimension 



19 



20 



2l|. 



Surprisingly, despite all of the work that has been performed on the phase diagrams with 
a strong-coupling analysis, there are only limited results for the momentum distribution 
functions. The first few terms of the structure factor have been determined to high order in 
„„e dimension □ a.d the .e.o d.t.but,o„ fu„ct,o„ has bee„ exa.i„ed ,„ o„e 



[l^ . A recent random phase approximation (RPA) has been carried 
which corresponds to the exact solution for the momentum distribution in the 



and two dimensions 
out 



infinite-dimensional limit (see also Ref. 



23l ). In this contribution, we present an alternative 



formulation of the strong-coupling perturbation theory for the many-body Green's functions, 
which can be immediately employed to evaluate the momentum distribution function as 
a power series in the hopping divided by the interaction strength for each value of the 



3 



momentum. Recently a similar strong-coupling formalism to ours has been proposed 2J] and 



used to calculate the momentum distribution in three dimensions through second order [25|. 
We take our strong-coupling expansion and, guided by the exact solution from the RPA, we 
construct an ansatz for the scaling behavior of the momentum distribution function and then 
employ it to produce analytic expressions for the momentum distribution that are accurate 
for all values of the hopping within the Mott phase. These results could prove useful as 
a simple means to check against experimental data on more recent Bose Hubbard model 
|. We also take the results for the scaling behavior of the momentum 



systems 



26 



27 



distribution and use it as a phenomenological ansatz for the scaling behavior of the phase 
diagram that sums many more terms than the original ansatz. Comparing that result with 
the QMC data in two and three dimensions also shows excellent agreement. 

We write the bosonic Hubbard Hamiltonian in the presence of a potential in the form, 



hap 
j 

Hoj = [Vrirj) - fi]nj + ^hj{nj - 1) (1) 

T-Chop = -^tjyajajv (2) 

id' 

Here j,j' label the sites of a (hypercubic) lattice in d dimensions, with a lattice constant 
which we set equal to 1 (the unit of distance); is the position vector of the j*'^ site 
as measured from the center of the system. The symbols and ajr denote creation and 
destruction operators for bosons at lattice site j. These operators obey the commutation 
relation [aj/,aj] = 6j'j] fij = a^^aj is the boson number operator at site j. VTivj) is the 
trap potential (which is usually assumed to be a simple harmonic oscillator potential) and 
the repulsive contact interaction is given by U . Note that the trapping potential could 
also represent a diagonal disorder potential, if desired, but we will not discuss that case 
further here. The chemical potential /i controls the average number of particles, ijji is the 
amplitude for bosons to hop from site j' to site j. We consider a general ijji for the formal 
developments we present in the earlier parts of the paper, but later specialize to the case of 
nearest-neighbor hopping only, with amplitude t, on a hypercubic lattice in d dimensions. 

As explained above, our aim in this paper is to calculate the properties of the Hamiltonian 
in Eq. ([2]), in particular, its momentum distribution function. The momentum distribution 
function is related to the atom-atom correlation function [see Eq. ([6])] involving atoms at 



sites j and j', which is given by 



(3) 



This expectation value can be calculated from the single-particle "thermal" or "Matsubara" 
Green's function, defined in the standard way 



as 



by choosing r = and r' = 0"^, the positive infinitesimal; i. e., 

Cj'j = -Gjj'{0,0^). 



(4) 



(5) 



Here, as usual, P = l/iksT) is the inverse temperature, < r, r' < /3 are "imaginary time" 
(henceforth ^H-time") variables, %- is the i-time-ordering operator, and Z = Tr[e~^^] is the 
partition function. The momentum distribution function measured in the time of flight 
experiments is proportional to the Fourier transform of the atom-atom correlation function: 



(6) 



where TV is the number of sites in the lattice (we do not discuss the proportionality factors, 
which arise from the Wannier wavefunctions of the trapped atoms, as that is not germane 
to the work we present here). 

Specializing to the case of nearest-neighbor hopping on a hypercubic lattice in d- 
dimensions, we report our main result which is the general strong-coupling expansion for 
the (T = 0) momentum distribution of the Mott phase with a density n up to third order 
in the hopping: 



n <{ 1 - 2(n + 1)^ + 3(n + l){2n + 1] 



uJ \u 



- 4(ra + l)[5n2 + 5n + 1] 



+ 



'-{n + l){26n^ + 26n + 5) 



-(n + l)(23n2 + 23n + 2) 



4d 



f/3 



f/3 



(7) 



where ek = — tX]^ 6xp[ik ■ 5] is the bandstructure, 5 is a nearest-neighbor translation vector, 
and d is the spatial dimension. 



Readers who are mainly interested in seeing how accurate this expansion is when apphed 
to exphcit cases, are encouraged to skip the next section which develops the formal techniques 
needed for obtaining the expansion, and proceed directly to Sec. Ill where we use the 
expansion to develop a scaling analysis and compare results to exact numerics. 

The manuscript is organized as follows: in Sec. II, we present the formalism for the 
strong-coupling expansion of the Green's functions and produce explicit formulas through 
third order for the one- dimensional lattice, the two-dimensional square lattice and the three- 
dimensional cubic lattice, along with the infinite-dimensional hypercubic lattice. In Sec. Ill, 
we present our scaling analysis for the momentum distribution in the first Mott lobe and 
compare those results to available numerical data from QMC and DMRG calculations; we 
also discuss the phenomenological approach to the phase diagram in two and three dimen- 
sions. Conclusions and a discussion of future directions follow in Sec. IV. Two appendices 
contain some of the more technical results. 



II. STRONG-COUPLING FORMALISM FOR THE GREEN'S FUNCTIONS 



The strong-coupling expansion we develop in this paper enables one to calculate Gjj' 
[in Eq. (jl])] and hence Cj'j [in Eq. ([5])] as an expansion in powers of Hhop [in Eq- (E])], 
with respect to regions of the system which are either normal or Mott-insulating (i. e., not 
superfluid ISw). For this purpose, we use the following standard relation 29| to define the 
i-time-ordered product for the evolution operator in the "interaction picture" : 



U{t, r')e 



U{t,t') = %exp[- J dTiHhopiri)] 



(8) 
(9) 



where, for any operator A, we define the time-dependent operator A{ti] 



on Wo 



Ae 



Using the properties of and the rules for composition for products of U, it is straightfor- 
ward to show that 



29| 



{X[UiP,0)a,{T)a\,{T')])no 



(10) 



The strong coupling expansion we use in this paper is obtained straightforwardly by expand- 
ing the exponentials in U [in Eq. (fTO!) ] in powers of Hhop and evaluating the resulting traces 
with respect to the equilibrium ensemble of Tio- The term of order m in such an expansion 
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for the numerator in Eq. (fTOj) is given by 

1 

x{T^[a,{T)a]jT^)a,,jT^) ■ ■ ■ a]^(r+)a,. (ri)a],(r')])Ho. (11) 



Since Tio, as defined in Eq. ([2]), is a sum of separate terms for each site, the thermal average 
in Eq. (fTTl) factorizes into a product of factors, one for each of the sites on the lattice, in 
terms of the multiparticle single-site Green's functions at these sites defined in the standard 



way 



as, 



^,(n,rO = -{%[a,{r{)a]{r[)])n,, (12) 



Q]\ti,T2]t'2,t[) = {Tr[aj{Ti)aj{T2)a]{T2)a]{r[)])no, (13) 



Note that these are total Green's functions, containing both connected and disconnected 
parts. Furthermore, each site that appears must occur an even number of times in the 
thermal average, half as indices of creation operators and half as indices of destruction op- 



erators [30|]. Similar considerations apply to the terms in the expansion for the denominator 
in Eq. (ITUl) [except for the absence of the operators aj(r) and a],(0)]. As we discuss in more 
detail below, the combination of the two expansions order by order leads to a cancelation 
of all "disconnected" terms, i.e., those involving products of thermal averages for clusters 
of sites that are not connected via hopping matrix elements to the sites j and j', as well as 
to the fact that the remaining terms can be written entirely in terms of the "connected" or 
"cumulant" multiparticle Greens functions, corresponding to the well known linked cluster 



theorem 



29| 



Using the above considerations, it is straightforward to write down systematically the 
terms in the strong-coupling expansion for Gjj'^T, t'). We denote the m*^ order contributions 
with a superscript (m). The different terms contributing in m*^ order can also be associated 
with "diagrams", which correspond to lattice "walks" or "world lines" for a particle which 
starts from site j' at i-time r' and reaches site j at i-time r after m steps (with each 'step' 
corresponding to a hop along the lattice, e. g., from site j[ to site ji induced by tj^jj, 
and in between the steps, the particle undergoes i-time "evolution" , which proceeds either 
forward or backward in i-time). These processes are shown graphically in Figs. [1] and [2l 



These diagrams are the strong-couphng analogs of the standard diagrams of many-body 



perturbation theory 29|], from which, after some practice, the terms can be written down 
by inspection. A p particle Green's function at a particular site appears when a walk visits 
that site p times. For m > 2, as we show below, the contributions can be classified further 
according to a hierarchy of decreasing powers of 1/z, where z is the coordination number of 
the lattice, by recombining contributions from intersecting and nonintersecting walks, and 
we denote these with further superscripts, as (m;0), (m; 1), etc. We give below the terms 
contributing to Gjj'{r,T') up to third order in Hhop, and their associated strong-coupling 
diagrams. 

The zeroth and first order terms are almost obvious. 

GfXT,r')=6,,g,{r,r'), (14) 

G'jjlir, r') = -t,,, / g,{T, n)gAn, r') ^ r'). (15) 

J Tl 

Here, and below, for notational convenience we denote integrals over i-times by integral 
symbols with subscripts, rather than by the standard notation. To second order, a 2-step 
lattice walk can either move to a distinct site two steps away or return to the starting site. 
Hence we get two terms from the numerator of Eq. ( fTOl) . the top equation is when the hop is 
to a different lattice site, the bottom equation is when the hop returns back to the original 
lattice site: 

Gf^^\T,T')num = {)-- 5jj')^ijnihj' I I (r, Ta) (xg, Ti) (n , t') , (16) 

- (17) 

G^f^ (r, T')num = Sjf tjj.tj^j / / (r, Tl, Tg, r')^ii (Ts, Tl) (18) 

jl JriJTi 

- (19) 

jl 

where the subscript num denotes that these are the terms coming from the numerator in 
the expansion for the Green's function. 

We have introduced a new notation above, letting ©j^j-/ denote the product of single- 
particle Green's functions at the sites that appear in the m-step lattice walk specified by its 
lattice indices, starting from right to left, together with the corresponding hopping ampli- 
tudes; the i-time arguments indicating the starting and ending i-time, the m intermediate 
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FIG. 1: (Color online.) Strong-coupling "diagrams" for the single-particle Green's functions up 
to second order in t. The horizontal directed dashed lines indicate the hopping matrix element t 
between the sites labeled, and the vertical lines indicate single-site Green's functions Q evolving 
between the respective i-times. The ellipses (yellow) at multiply visited sites denote the appearance 
of connected or cumulant n-particle Green's functions. 

i-times being integrated over, ^j-"!^/ is defined similarly, except that it necessarily involves 
self-intersecting lattice walks where one or more sites are visited multiple times, and the 
product involves r-particle Green's functions at a site that is visited r times, with the inter- 
mediate i-time arguments being determined by the sequence specified in the lattice walk. For 
any given m-step lattice walk both ©("^^ and can clearly be written down by inspection. 
To correctly obtain obtain G-L we need to subtract from the above two terms the term 
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FIG. 2: (Color online.) Strong-coupling "diagrams" for the single-particle Green's functions cor- 
responding to third order in t. The horizontal directed dashed lines indicate the hopping matrix 
element t between the sites labeled, and the vertical lines indicate single-site Green's functions Q 
evolving between the respective i-time. The ellipses (yellow) at multiply visited sites denote the 
appearance of connected or cumulant n-particle Green's functions. 

that arises as the product of the second-order contribution from the denominator of Eq. (fTOj) . 
corresponding to closed loop lattice walks involving the sites j and ji, given by 



^(2) ^ , 



J TO. J T^ 



(20) 



and the zeroth-order term from the numerator, namely G^^I{t,t'). The net result for G*-^-* 
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can be reexpressed as the sum of the following two contributions: 

G^iry) = E^Su' (21) 

ji 

Jt'2 J Tl 

^^n'Y.^.y-') (22) 

h 

= ^n' E{^S.(^' ^') - ^S[^y - ^(^' ^')^S}. (23) 



n 

Here 

(r, n; r2, r') = [^f (r, n; r2, r') - G.ir, T2)g,in, r') - Q.ir, r')^,(ri, Ts)] (24) 

is the "cumulant" or "connected" part of the two-particle Green's function at site j. 0^"^^ 
is defined similarly to ^^"^^ except that the multi-particle Green's functions that appear in 
^j,g QQnnected Greens functions. Note that the prefactor ^ present in Eq. (fTTj) no 
longer appears in the above equations, as it has been canceled by the 2! ways of choosing 
the two distinct hopping matrix elements in the expansion. 

Similarly, the third-order contributions involve three-step walks. From the numerator of 
Eq. fllOp we get the following terms according to the types of walks involved. 

G^jj' )num = —(1 — Sjj') ^^(1 — <^i2j')(l ~ ^jji) ^jj2^j2jl^jij' 

j2,jl 

X 



/ / / gAr,r3)g,,{T3,T,)g,Ar2,n)gAri,T'), (25) 

J T3 J T2 J Tl 

^ii'^^(^' ^')n«m = -(1 - Sjj') ^(1 - Sjji) ^jj'^j'jAjij' 

ji 

X / / / ^j(r,r3)^j,(r2,ri)^j/(ri,r3;r2,r') 

-(1 - 6jy) ^(1 - 6j^j>) tjj^ij^jtjf 
h 

X / / / gY{T,T2,T3.n)g,,{T3,T2)g.M^T'), (26) 

J T?. J T2 J T^ 



' T3 J T2 J Tl 

.(3;c) , ,s 1 



Gjj' i'^j'^)num — 2] (l ^ii') tjj'tj'jtjj' 



X 



/ / / 0i\T,T2;Ts,n)g'/{rs,TuT2,T'). (27) 

J T'i J T7 J T^ 



' T3 ^ T2 -J Tl 

.{3;c) 



Again, in all cases except for the case of Gnum^ the 3! ways of choosing the three distinct 
hopping matrix elements involved completely cancels the ^ in the expansion. In case of 
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Gnum, two of the hopping matrix elements are identical, so they can be chosen in only || 
ways, hence there is a factor of ^ left uncanceled. 

The restriction (1 — 5jji) in the third-order contributions above is redundant except on 
nonbipartite lattices, such as nearest-neighbor hopping on a triangular lattice, or on a hy- 
percubic lattice with second-neighbor hopping, where one can return to the starting site 
after three hops. In such cases, one has the additional term 



X / / / 0'/ir,n;n,T')gj,in,T2)g,,ir2,n). (28) 

Note that in this term the constraints on ji and j2 are actually redundant and can be 
omitted. 

As in the second-order case, the above third-order terms can be recombined with the 
terms that arise as products of the first-order term from the numerator and the appropriate 
second-order terms from the denominator in Eq. (fTOl) and reexpressed compactly in terms 
of the connected Green's functions (3^^^ and (3^^\ One gets 



G^i-y) = -E^gLu'(-'-'), (29) 

12, jl 

ii j2 
Gp^{r,r') = -^^%,{ry). (31) 

In nonbipartite cases, one has to add to this the additional contribution 

G^f (r,r') = -5,,,E^SU.(-'-'). (32) 

32 ,31 

Note that in all the cases, use of the connected Green's functions allows one to avoid the 
clumsy restrictions on the intermediate sites that need to be summed over, and in addi- 
tion, automatically includes the terms contributed by the denominator of Eq. ffTOj) . The 
diagrams that represent the above are shown in Figs. [T] and [2J The same results can also 
be derived using more formal methods, such as functional integrals, generating functionals, 
and functional derivatives, but we do not go into such details here. 

Next, we discuss the evaluation of the multiparticle single-site Green's functions at a site 
i as defined in Eqs. f|T2|) and f[T3|) . The eigenstates of Tioj in Eq. ([2]) are also eigenstates 

12 



of the number operator nj, and can hence be labeled by positive integers n = 0, 1,--- 
corresponding to the number of bosons at site j, with energy eigenvalues which we label as 
ej^n- One has, 

T<oj\j,n) = ej>|j,n); ej> = [Vrirj) - ij]n + ^n{n - 1). (33) 

The partition function of the j^^ site, and the Boltzmann probability of occupancy of |j, n) 
in the thermal ensemble corresponding to Ti-Qj, are given respectively by 

Zj = ^exp (-/5ej-„); pj-„, = exp (-/3ej-„)/Zj. (34) 

n 

It is convenient to define the ladder operators 

X+^=\j,n + l){j,n\, Xr^ = \j,n-l){j,n\. (35) 

One can easily see that 

(^) = E = E ^^n, (36) 

n n 

where 

are the "particle" and "hole" "excitation energies" (with respect to the state with n bosons 
at site j) induced by the ladder operators and '^j~„j., respectively. 

Using the above, and the rather obvious rules for products of the ladder operators, it is 
easy to verify that, for the 1-particle Green's function, we have 

Gjin, r2) = E PjAi^ + 1) e^"^-"^)^-" e{n - r^) + n e'-^^-^'K- Bin - n)]. (38) 

n 

There is a compact way of working and writing this out which easily generalizes to n- 
particle Green's functions. Let P label the 2! possible permutations of (1,2), corresponding 
to (1,2) — > (P1,P2), and ^j(P) denote Qj{Ti,T2) in the domain (rpi > rp2)- Furthermore, 
define cti = — 1 , (T2 = +1 ; tf^ = e^„; and X^J;; = X^^. Then, one has. 



n ni,n2 



=1,2 



E ^i'" L^£l Jn + exp [rp2ej;^ - rpie^- ^" ] . (39) 
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As is easily verified, for tlie identity permutation, corresponding to PI = 1,P2 = 2, tliis 
reproduces tlie first term in Eq. fl38l) : for tlie permutation corresponding to PI = 2,P2 = 1, 
it reproduces the second term in Eq. (155]) . 

Now, for tlie case of the 2-particle Green's functions, let P label the 4! possible per- 
mutations of (1,2,3,4), corresponding to (1,2,3,4) (PI, P2, P3, P4), and g^/(F) denote 
Qj^Ti, T2;ts, T4) in the domain (rpi > rp2 > rp3 > Tm)- For this case, we define ai = (T2 = —I 
and (J3 = (74 = +1. Using these definitions, we can show that 



ni,--- ,n4 



exp rp£e^. „^ W + — ^ — A^- „^ 

i=l,---,4 



n 



El 1 - an 1 - o-p2 / , 1 + o'ps / 1 
Pj,n \ \ n- cTpi H \ n + H \ n + 



2V 2V 2V 2 



X exp [rp46^- + rpse--^^^^) - rp^e-r^,, - rpi67-«]. (40) 

For example, in the domain (ri > T2 > T3 > T4), corresponding to the identity permutation, 
this formula gives, 

= {aj{ri)aj{T2)a]{T3)a]{T4))no, 

= Yl + + 2) exp [{n - ri)e+„ + (rg - T2)e+^^^]. (41) 

Using the above results, we can now readily compute the terms in the strong-coupling 
expansion of Cj'j in Eq. ([5]) up to third order in the hopping amplitude. In the equations 
below, we denote (E^!], ^ ^^'"^^(O, 0+), = ©5:^^,(0,0+), and ^^'"J, = ^^'"^^(O, 0+). 

The zeroth order term is 

Cf:] = -G^iO, 0+) = 6,,y{n,)n,^ = 5,,, ^ n p,,„. (42) 

n 

The first order term is 

^r] = 0^) = / ^.(0' n)0An, o+) = <E^. (43) 

J n 
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The i-time integral is straightforward to evaluate. Using Eq. (138!) or Eq. (139|) . we find, 

1 - exp [-/9(eT^ + e+ „,)] 



f ~ -I- 

j,n ' j',n' 



^jj'^^i^' + 1) Pj,nPj',n- 
n,n' 

i^^, J2 n{n' + 1) I -^^'^ ■ ^^'"-^ 



t • ^p- I ^(^' + ^) + (n + l)?!' I ^^^^ 



This can be represented by the diagram labeled C^^^ in Fig. [31 The third line of Eq. 
is written in a form that can be directly constructed from this diagram in a way that is 
immediately generalizable to higher order (see below). The fourth line contains a second 
form of the same result, obtained by relabeling the bosonic occupation numbers in the second 
term of the third line in a way that makes the zero temperature limit obvious. 

There are two second order terms in Cj/j corresponding to the two terms in Gjjr [Eqs. ( !2T|) 
and (|23|)]. 

cjf = -G<f(0,0-) = -^4|,. (45) 



Jl 



and 



jl 

= -^n' - - (46) 

jl 

Similarly, one obtains the various terms contributing to Cj'j in third order, which we label 
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Cpf \ Cff, ■ ■ ■ , by setting r = and r' = 0+ in Eqs. ((29]) , ([30D , ■ • ■ . One finds 
^^(3,0) _ sr^A^) 



j'j jj'jij' Z— ihjj' 

jl 32 



/L^L n ji] 33 nf 33 3 31^ 



31 



^(3,2) ^ J_^{3) 

j'i 2! 

= . , _ 201 . J,- - 20^ . .,l„-,v + 2C^.^.) . ., + 2€.^^lz^^h. (49) 

2\'- 33 33 1- 33 313 '^^-^ ^ 33233^^^-^ ' Jj'jj' ' 33' 3' 3 i ^ > 

The i-time integrals that appear in these expressions are most conveniently evaluated by 
splitting them up into separate integrals corresponding to each of the different (m!) i-time 
orderings of i-time integration variables (in m*'^ order). With each such i-time ordered term 
one can associate a unique diagram, as shown in Figs. [3H6]up to third order. The diagrams 
are labeled by the sites that appear, the "initial" (= "final") and the "intermediate states" at 
these sites as determined by the boson occupation numbers at these sites in each of the i-time 
intervals, whose labeling corresponds to the boson creation and destruction processes at the 
sites. (The boson occupation numbers at the sites that do not appear in a diagram do not 
change with i-time, and play a spectator role, and hence do not appear in the contributions 
to Cj/j.) The "matrix elements" that are associated with these processes are then uniquely 
determined and can be written down by inspection from the labeling. For such a diagram 
of m}^ order, let and , ■ ■ ■ , denote the energy eigenvalues of for the initial 
(or final) state and the m intermediate states respectively. Then the i-time integral is of the 
form, 

^ Jo Jo Jo 

This is easily evaluated using Laplace-transform techniques, as shown in Appendix [X] If the 
energies are all distinct, then one finds that the integral is the following sum of m + 1 terms. 

ImiP'i ^-Qm) ■ ■ ■ ) £ai, £ao) = / ^ J_ J_ TT; T- (51) 

i=0 I'i^i ^ "^l' 
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Note that only energy differences appear in the energy denominators in this expression, and 
they are related in a simple way to the boson creation and destruction processes at the sites 
that appear in the diagrams; these can be written down by inspection from the labeling 
shown in each diagram. As the initial and intermediate states at all the sites that do not 
appear in the diagrams are constrained to be the same, one can replace the Boltzmann 
factors for the initial and intermediate states by a product of the density matrices for just 
the sites that appear in the diagrams. 
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FIG. 3: Strong coupling "diagrams" for the correlation functions up to second order in t. The 
horizontal directed dashed lines indicate t between the sites labeled, and the vertical lines indicate 
single site Greens functions Q. The ellipses (red on line) at multiply visited sites denote the 
appearance of connected or cumulant Greens functions. Initial and intermediate state labels for 
the different possible i-time orderings shown are also indicated. 

The expression in Eq. (15T!) is non- singular and remains well defined even when one or 
more of the energies Sao,Sai, ■ ■ ■ , £am become equal, as clearly happens, for example, in the 
diagrams for C^"^'^^ (see Fig. [3]). For example, if one and only one pair of energies are equal, 
say, = Sap , then instead of Eq. fl^ one should use the expression 



E— Il ig _g ) +— I/^-Ett^TI' n (g .£ y (52) 
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FIG. 4: (Color online.) Strong-coupling "diagrams" for the correlation functions C^'^'^\ The 
horizontal directed dashed lines indicate the hopping matrix t between the sites labeled, and the 
vertical lines indicate single-site Green's functions Q. The ellipses (red) at multiply visited sites 
denote the appearance of connected or cumulant Green's functions. Initial and intermediate state 
labels for the different possible i-time orderings shown are also indicated. 

which it reduces to in this limit (see Appendix |A] for further details). 

The diagrams labeled C^^'°^ in Fig. [3] shows the two diagrams for corresponding to 

the two i-time orderings T2 > ti and ti > T2 ■ Each diagram gives rise to three contributions 
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FIG. 5: (Color online.) Strong-coupling "diagrams" for the correlation functions C*-^'^^. The 
horizontal directed dashed lines indicate the hopping matrix t between the sites labeled, and the 
vertical lines indicate single-site Green's functions Q. The ellipses (red) at multiply visited sites 
denote the appearance of connected or cumulant Green's functions. Initial and intermediate state 
labels for the different possible i-time orderings shown are also indicated. 
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FIG. 6: (Color online.) Strong-coupling "diagrams" for the correlation functions C*-^'^^. The 
horizontal directed dashed lines indicate the hopping matrix element t between the sites labeled, 
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sites denote the appearance of connected or cumulant Green's functions. Initial and intermediate 
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as pointed out above, which are easily written down using the above rules, leading to 



n,ni,n' J'" ''ji,niJ\'^j,n "T "^/.tiV 

I Pj,n-1 Pji,ni Pj',n'+1 Pj,n-1 Pji,ni+1 Pj' ,n' 



Pj,n-1 Pji,ni Pj',n'+1 _| Pj,n Pji,ni-1 Pj',n'+1 j | ^^g^ 



+ ^ji,ni)i^t,n-l ~^ ^j',n'+l) i^j,n + ^ji,ni-l) (^ji,ni-l + ^■',n'+l) 

^ n(ni + l)(n' + l) 

'-jjihij' / J Pj,n Pjx,nx Pj' ,n' \ [ 



(n + l)(ni + (n + l)ni(n' + 1) 



^ |. nriiiji' + 1) 

(n + 1)^172^ n{ni + l)n' 

Again, the second form of the result, Eq. flM|) . is obtained by appropriately relabeling the 
bosonic occupation numbers in four of the six terms in the first form [Eq. fl53|) ]. and is easier 
to use at T = 0. 

Similarly, from the diagrams contributing to cgj^- (labeled C(2.i)) and zgj (labeled Z^^)) 
shown in Fig. [31 and using Eq. we obtain, for the two equivalent forms for each. 



n.ni 



+ (n + Dm |-^i^^^(/3 - , . ' - J + fa.^..- ] } (55) 



ni X 



n.ni 



^ n{ni + 1 ^ ^ {n + l)ni{n + 1)^ 



J," ' n,niJ 



+ +6T " (e+ +67 )^ + (eT +e+ )2 ^ ^ ^^^^ 
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+ f 4- f"^ -I- f ~ P= 

n,ni ^"J'" ii,"!'' i," ' ii,"!'' ^ 3,n-l ' ^ji,ni+l/ 

+ (n + [ 'r'^p - T-T^—^) + J^""i^r-\. ] } (57) 

V j,n "T ^ji,n-i) y j,n ' ji,ni) \^j,n+l ' *^ji,ni-l) 



^jji ^h 'j Pi-.n Pji ,ni 



X 

n,ni 



n(ni + l) 1 ^ + 

Simplifying these expressions we find 

+ )2 - )J > (59) 

-7(2) t i J /5 r ^(^1 + 1) , (^ + 1)"-! 1 X /^n^ 

^jn = ^ni^hj 2^ Pj,n Pjum { P [ . + ^ + J }• (60) 

n,ni yj,n~^ ji,ni) \'^j,n~^^ji,ni) 

Finally, Figs. HI [5] and [6] show the diagrams corresponding respectively to ^fl^ ,v, ^f}^ 

— /o\ 

(the digrams for ^^jj^jji can be obtained from those in Fig. [5] by symmetry and a simple 
relabeling) and ^jj/jj/- (For simplicity, since we do not discuss nonbipartite lattices in detail 
in this paper, the diagrams for ^fj^jij shown.) In each case, there are 6 possible 

orderings of the i-time variables t^, T2 and Ti [listed in the order (ts > T2 > Ti), > t-^ > Ti), 
{t2 > Ti > T2), {r-i > Ti > T2), {ti > T3 > T2), and (ri > r2 > Ts) below]; and from each i-time 
ordering we get four contributions corresponding to the initial and three intermediate states 
(apart from the subtractions arising from the connected two-particle Green's functions). 
The contributions can be written down straightforwardly using the rules stated above, and 
we present them in Appendix [Bl 

It is easy to see that the methods we have discussed above permit one, in principle, to 
similarly write down the contributions to Gjj> and Cj>j to higher orders as well, though the 
calculations will become increasingly tedious unless one can find a way to automate them. 
However, it is possible to calculate sums of subsets of these contributions to arbitrary orders. 

The easiest subset to sum is G^J^,'^\ with contributions coming entirely from "self-avoiding 
lattice walks" while computing thermal averages, and therefore involving only single-particle 
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single-site Green's functions, but ignoring the self- avoidance constraint while summing over 
the different possible walks. From the above analysis it is clear that the resulting term is 

--<(m;0)/ /N _ i i , I ! ■■■ ! I 

T,n J Tm-1 J T2 J T\ 

X Gj{r, Tm)Qj^_^{Tm, Tm~l) ' " " Qh{T2, t')- (61) 

The sum of these to arbitrary order, together with G^J] [Eq. (fTSl) ] and G^y [Eq. (fT4|) ]. 
correspond to a geometric series for the Green's function regarded as a matrix (denoted by 
bold- face letters and/or square brackets below) with lattice sites and i-times as indices, and 
corresponds to the well known "Random Phase Approximation" (RPA) result, 

[G''''%r{r, r') = [Gjr\r, r')6,y + 5(r, r%y. (62) 

By taking a Fourier transform with respect to the even Matsubara frequencies iVtm = 
2m7TkBT, m = 0, ±1, ±2, ■ ■ ■ , one can write this as a matrix equation involving only lattice 
indices: 

[G«^^]7;(zi]™) = [g,{zQ^)]~%,, + t,r, (63) 

where, from Eq. fl38l) . the single-site Green's function at a fixed frequency is easily obtained 

as 

e,(,:n„.) = i^J^ - j^^^l (64) 

The RPA correlation function can then be straightforwardly obtained, in view of Eq. ([5]), as 
the Matsubara frequency sum 

Cf, = -[G],y (0, 0+) = -J2 G.^itQ^e'"'-''' , (65) 
by using [G'^^^] for [G]. While the frequency sum can in principle be evaluated using 



standard contour integral techniques 29|, in the inhomogeneous case the above calculation 
involves a matrix inversion with respect to the site indices. 

The calculations simplify, however, for the case of a homogeneous system, i. e., without 
a trap or disorder potential, and for T ^ 0. Then all the sites are identical, and in a generic 
case, the ground state of TYqj at every site corresponds to the same fixed boson occupancy 
which we denote n. In this limit, one has e„ = — /i?T,+ ^n{n — 1) whence = — // + Un and 
e~ = yU — IJin — 1). The RPA momentum distribution can be calculated exactly analytically 
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in this discussed in detail by Sengupta and Dupuis [2^. From Eq. (16^ . at T <^ 

we get, for all sites j, 

(n + 1) n 



[iVt^ - e+) {iVt^ + e-; 



(66) 



Hence, 



1 - , Z]^ 



(67) 



with poles at = [e^ + — e~ ± -^/e^ + 2e]JJ{2n + 1) + U'^]/2 and residues determined 
in terms of = (-E,^ + /i + f/)/ (-E,^ — -E^ )• The RPA momentum distribution at T = is 
just the negative of the spectral weight of the pole at E^; i. e., 

= ^k - 1 = , K + l^ + U (68) 

The challenge, of course, is to go beyond the RPA. One way to achieve this, by summing 
further infinite subsets of contributions to Gjji beyond the RPA, is by using the Dyson 
equation (compare Eq. [62]) . 



[G]^'(^, r') = m-\r, r')b,,, + 5(r, r')i,,, - sg)(r, r') - Sg)(r, r') - • • • , (69) 

where S*^'") denotes a self-energy correction that corrects RPA to order t™. By re-expanding 
the inverse of this equation and comparing with the expansion for G^ji discussed earlier, it 
is straightforward to obtain the following expressions for the self-energy corrections up to 
third order in t. 



/ / Kl-'(r,T2)GP;"(7.,,T0|e/|-'(Ti,r') 

J T-2 J Tl 

j j E[^^r'(^'^2)^g,(r2,ri)[^,r^(ri,r'), (70) 

J TO J 



with 0^^) as given in Eq. fl2^ . Similarly, 



Jt2 Jtx 

J Tl J TO 



24 



It is straightforward to verify using the expressions given in Eqs. ( !23|) and (l30l) - (l32l) . that the 
term involving G -] exactly cancels the two terms involving G -', above, and one obtains, 




(71) 



One can in principle evaluate these expressions for the self-energies explicitly using the tech- 
niques discussed above, and thereby determine spectral functions as well as the momentum 
distribution function using Eq. fIBSl) . We plan to complete such work in the future. 

However, in this paper we adopt a different procedure for calculating the momentum 
distribution function for the homogeneous case and in the T ^ limit. We directly evaluate 
the expressions for Cjij up to third order in t. Then we use a scaling ansatz for the momentum 
distribution function determined in such a way that when expanded in powers of t, it agrees 
with our calculated results, thereby effecting an infinite order resummation in a different 
way which automatically has the correct critical behavior at the Mott-superfluid transition. 

The direct evaluation of our expressions for Gjij in the homogeneous, T ^ limit is 
straightforward. The T — > limit is easiest to implement using the second form of these 
expressions, where, just as discussed above in case of the RPA, for T <^ the sums 
over the initial and intermediate states are all restricted to n. The excitation energies 
that occur in the energy denominators in these expressions are given by (e^ + ^n) — ^ i 



find 



^n — 1 



■ e+ + e" 

n ' n 



3U and e" 



'n — 1 



4U. Hence, we 



2n{n + 1) 
U 



(72) 
(73) 



For the second-order terms we obtain 



i2n 



r{2) 
"jjij 



?>n{n + 1 
3n(ri + l)(2n + 1 



3n2(ri + 1) 
IP 



[/2 



2n(n + 1) 
U 



(74) 
(75) 
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Hence, using Eqs. ( H5!) and (H6|) . we find 



(2,0) 



3n(n + l)(2n + 1) 
t72 



(2,1) 



A r^(2,0) 



(76) 
(77) 



Note that in C^"^'^^ the (divergent) temperature dependent terms from C'-^-' and Z^"^^ exactly 
cancel, as they ought to. 

Finally, we consider the various third-order terms. We get the following results 

"4n(n + l)^l \3n'^{n + lf] \4n^{n + iy 



r{3) 



t/3 



+ 4 



f/3 



t/3 



4n(n + l)(5n2 + Sri + 1) 



f/3 



whence, 



Furthermore, 



^(3) 



(3,0) 



4n(?i + l)(5n2 + 5n + 1) 
t73 



jj'jlj' j^' j'-jl Jl-?' 



2n(n + 1)3 ^ n^(n + l)^ ^^^^ _ 



f/3 



f/3 



2n2(ra + 1)^ _^ ra(ra + l)(n2 - 1) ^ 72(72 + l)(n^ + 2n) 



f/3 

722(72 + 1)2 
f73 



{PU - 2) + 



3f/3 3f/3 

72(72 + 1)(722 — 1) 7i(72 + l)(7l2 + 272) 



3t/3 



+ 



272^(72 + 1) ^ n\n + lf ^^^^ _ 



f/3 



[/3 



222(22 + 1)(4722 + 4r2 + 1) ^ ^722(72 + 1)^ ^ 



3t/3 



Similarly, we find 



r(3) 



272(72 + 1)(4722 + 472 + 1) ^ 72^(72 + 1)2 



3f/3 



t/2 



Using these and Eq. (148|1 . we get 



33 



■33 



31 



32 



272(72 + 1)(26722 + 2672 + 5) 

31P 



(78) 
(79) 

(80) 



(81) 
(82) 



(83) 



(84) 
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Next, 



^(3) 



2 ^ ..T ^ {pU - 2) 



+ 4 



n {n + 1 



f/3 

-(/5t/-2) 



t/2 



/3 



4(n - l)n(n + l)(n + 2) 
2ra(ra + l)(7r22 + 7^ + 2) 

[73 



Hence, using Eq. ( 149|) . we obtain 



(-7(3>2) 



n(n + l)(23n2 + 23ra + 2) 
3f73 



(85) 
(86) 

(87) 



Note, again, the exact cancelation of the divergent temperature-dependent terms above. It 
is straightforward to verify that the Fourier transforms of the expressions for above 
up to third order agree with those obtainable by expanding the RPA expression in powers 
of ek (see below). 

III. SCALING ANALYSIS 



In the rest of this manuscript, we specialize to the case of nearest-neighbor hopping on 
a hypercubic lattice in d-dimensions. Combining the different contributions for Cj>j for 
such a lattice and Fourier transforming to momentum space, we arrive at the starting point 
for the scaling analysis, which is the strong-coupling expansion for the zero-temperature 
momentum distribution truncated to third order in the hopping and shown in Eq. ([7]). It 
is more convenient to reexpress the results for different cases in terms of the dimensionless 
parameters x = di/U and = ek/2(it. If we further consider only the n = 1 Mott insulator, 
we find 









-32 



22^^ 



19 



2 

rf2 



X 



I. e., 



nk = 1 - S^ks; + T^iW - 704^^x 
in infinite dimensions where x remains finite as — oo. 



3^3 



55, 



nk = 1 - sax + 12[6e^ - l]x2 - 32[22e^ - -]i^.x\ 



in three dimensions, 



nk = 1 - S^kX + 18[4^^ - l]x^ - 32[22e^ - 9\i^x\ 
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(88) 
(89) 

(90) 
(91) 



in two dimensions, and 



nic = 1 - 8^2: + 36[2^^ - l\x' - 32[22^^ - l?]^^^ (92) 

in one dimension. Note that because integrals of odd powers of over momentum vanish, 
we can use the fact that the integral of the square of .^k over k is equal to l/{2d) to show 
that the integral of the strong-coupling expansion for nk over k is always equal to 1, as it 
must be. 

We show comparison of these truncated third-order strong-coupling expansions directly 
with exact numerical results and other analytic approximations below. It turns out that the 
truncated strong-coupling expansion does not work so well for the momentum distribution 
once the hopping is on the order of one fourth of the critical hopping for the Mott to super- 
fluid transition in two and three dimensions (and is even worse in one dimension). Hence, 
we use additional knowledge about the momentum distribution and how it scales near the 
critical point, along with the exact solution in large dimensions to create a phenomenologi- 
cal ansatz for the momentum distribution which produces analytical expressions useful for 
direct comparison with experiment. 

We start our scaling analysis with a general discussion. The zero momentum distribution 
function becomes critical at the critical value of the hopping for the Mott insulator to 
superfluid transition (called Xc). The critical behavior goes like nk=o — > ^(^~'') where ^ oc 
l/{xc — xY is the correlation length of a ci -|- 1 dimensional XY model [2! and r] and u are 



critical exponents in the usual notation [3l|. In two and higher dimensions, the correlation 
length diverges as a power law. The critical exponents for the two-dimensional Bose Hubbard 
model, which correspond to the three-dimensional XY model, are 77 = 0.04 and v = 0.67, so 
(1— //)// = 7s = 0.64. In three and higher dimensions for the Bose Hubbard model, the critical 
exponents are mean-field like, with 77 = and u = 0.5, so (1 — rj)]/ = 7s = 0.5. The one- 



dimensional case has Kosterlitz-Thouless behavior [3^, where 77 = 0.25, and the divergence 
of the correlation length has a Kosterlitz-Thouless exponential form C, oc ex])[W/ ^/x^^^^lE\, 
with Xc the critical point for the Mott insulator to superfluid transition. 

This critical scaling behavior does not provide enough information for us to determine an 
ansatz for the momentum distribution function over all momentum, because the distribution 
function is not critical for nonzero momentum. We use the exact solution in the infinite- 
dimensional limit, as given by the RPA solution, to guide us in how to proceed to develop 
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an appropriate scaling ansatz. The RPA form for the momentum distribution function, as 



22| 



discussed above, and reexpressed in terms of .^k and x, is given by 

2 v^l + 4(2n + l)aa: + 4e^x2' 

and this is the exact solution in infinite dimensions. A quick examination of the strong- 
coupling expansion for arbitrary dimensions, shows that the 0(1) terms are the same for all 
dimensions, when expressed in terms of x and ^k, and it is only the l/c?" corrections that 
differ for the different dimensions. Hence, the power-series expansion of the RPA form must 
produce all of the 0(1) terms. In finite dimensions, only l/rf" corrections are allowed. This 
motivates the following scaling ansatz for the momentum distribution function in two or 
higher dimensions (on a bipartite lattice) 

2 [1 + 2aekX + 46e^a;2 + + 8Je^x3 + 2f ek2;3]7« ' ' 
with d the spatial dimension. Note that in three dimensions, since 7s = 0.5 which is the 
same power law as in infinite dimensions, we must have d = 0. We will see this occur in the 
analysis below. 

In order to determine the parameters in the scaling ansatz, we propose three requirements 
of the formula in Eq. 0941) : (i) the power-series expansion of the scaling ansatz, in powers of x, 
must reproduce the strong-coupling expansion through the given order (in our case through 
third order) as shown in Eq. ([7]); (ii) we choose the scaling form to have the exact critical 
point Xci as determined by QMC, DMRG, or scaling results of a strong-coupling expansion 
for the phase diagram; and (iii) we require the integral of rik over all momentum to give 
n, the density of the bosons in the Mott phase. In two and higher dimensions, these three 
requirements will determine all of the parameters, which we now show; in one- dimension, 
we use additional information to determine the Kosterlitz-Thouless constant W , which then 
allows us to determine the complete scaling form. 

We begin with the infinite-dimensional case where the hopping scales like l/rf so that 
X is finite, but the coefficients c, c', e, and e' all vanish because they are 1/d? corrections. 
Expanding the scaling ansatz in a power series in x and equating the coefficients of the powers 
of X with the strong-coupling expansion in Eq. ([7]), yields the following: a = (2n + l)/'^s = 
2(2n + l),b = l/27s = 1, and d = 0, so we recover the RPA result in Eq. (p3|) . Since the 
critical behavior occurs at the point where 1 -|- 4(2n -|- l)C,kXc + 4:^]ixl = 0, and we evaluate 
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for k = 0, where .^o = —1, we immediately find that Xc = (n + 1/2) — y^n{n + 1), which is 
the exact critical point 0, 9] for all n. Hence, one can see that this approach automatically 
produces the right behavior for the large-dimensional limit. 

Note that the curvature of the RPA momentum distribution, with respect to is always 
one sign. In the truncated third-order strong-coupling expansion, the curvature of the 
momentum distribution function changes sign at = 1 when x ~ 0.034. This effect occurs 
for all finite dimensions as well. The scaled results, that are shown below, do not have a 
change in the sign of the curvature, and we expect that this does not occur in any of the 
exact solutions of the Bose Hubbard model. 

Since the momentum distribution function depends on the correlation length at k = 0, 
and so does the phase diagram, it is interesting to try a phenomenological exercise, where we 
take the critical behavior determined via our scaling approach for the momentum distribution 
and relate it to a determination of the phase diagram. Since we have a power law of a 



polynomial, instead of the simplest scaling dependence, which would go like -v/x^— such 
an approach is similar to summing an infinite number of terms in the expansion for the Mott 
phase lobes in the phase diagram. As an example, we make the following scaling ansatz for 
the Mott lobes 



/i 

U 



= n + A{x) ± [scaling polynomial] (95) 



where the scaling polynomial is the polynomial used for the momentum distribution at k = 
[which is 1 — 4(2n -|- l)x + for the infinite-dimensional case], and A{x) and B{x) are 
polynomials in x. Fitting the parameters to the third-order expansion for the Mott phase 
lobes, we find for the infinite-dimensional case that 



U 

which is the exact solution 



n---x±- v^l -4(2ra + l)x + Ax'^, (96) 



In finite dimensions, we will consider only the n = 1 case, because good numerical data 
is available for the Mott phase boundary in one, two, and three dimensions and we want to 
ensure that we produce the correct critical point Xc- We start with the three-dimensional 
case. Taking n = 1, we expand Eq. ( 194|) in a power series in x, and compare with the 
strong-coupling expansion in Eq. (l90l) . We find a = 6, 6 = 1, c = 144 -|- 4c'/3, = 0, and 
e = 224/3 + 58c'/9 + 4e'/3. At this point, the constants c' and e' are not determined. We 
fix e', by requiring the momentum distribution at k = to diverge at the critical point 

30 




FIG. 7: (Color online) Momentum distribution function for the three-dimensional case with x = 
0.0625 as a function of the band energy ek- Note how the QMC data agrees better with the scaling 
theory results than it does with the strong-coupling results or the scaled RPA, although deviations 
can be seen in the data. 



Xc = 0.10224 as determined by QMC simulation ITj. This produces the equation 



1 - 12a:, + ( 20 + -c'] xl - ( — + —c' + -e' ] xl = 0. (97) 
V 27 y ^ V 27 81 27 J " ^ ^ 

Setting Xc = 0.10224, and solving for e', yields 

e' = -122.2743 -f- 0.0571205c'. (98) 

(Note that if we instead set the coefficients c' and e' to zero, then the critical point would lie 
at Xc = 0.09805, which is about a 4.3% error.) The coefficient c' is determined by requiring 
J d?kn\^ = 1; we find that c' ranges from at x = out to d = —1.86 as x — * Xc- A simple 
polynomial fit to the behavior of c'{x) is 

c'{x) = 0.017166 - 0.71982a; - 161.093a;2 - 109.614x1 (99) 

We compare the strong-coupling perturbation theory to numerically exact results per- 
formed with world-line quantum Monte Carlo simulations of the Bose Hubbard model that 



employ the directed- loop algorithm 33|], in particular, its continuous- imaginary-time vari- 



ant [3J] . We have further improved the algorithm by omitting one-site vertices corresponding 
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to the U-term 



351 ] and also two-body vertices corresponding to the hopping term [36|. The 



latter modification is useful in reducing the memory and was crucial in the present simu- 
lation of the largest system {L = 64). The accuracy of the method is tested by comparing 
with exact diagonalization for small systems, and verifies the critical exponents with known 
results for the d+ 1-dimensional XY model. To further test that the true equilibrium distri- 
bution is sampled on large systems, several independent runs with varying lattice sizes are 
carried out, showing no systematic deviation, thereby ensuring that our numerical results 
are "exact" except for statistical errors. This QMC approach has already been applied to the 
problem of determining how the momentum distribution changes when the system becomes 
superfiuid fisl . 

We also compare the momentum distribution to RPA results. Since the RPA has a 
critical value of x that is smaller than the true critical value in finite dimensions, we plot the 
RPA results in Eq. (!93|) at a rescaled hopping value, corresponding to the same fractional 
amount of Xc- Namely, we choose xrpa = 0.0857864x/a;c((i). We call this the scaled RPA 
momentum distribution. 

The scaled results of the strong-coupling perturbation theory fit the numerical QMC data 
quite well. We compare with data at x = 0.0625 and x = 0.09 in Figs. [7]and[8|, respectively. 
The QMC data is for a 48 x 48 x 48 lattice at a temperature T = O.lt (T = 0.025t for 
X = 0.09); in all cases, we have carefully checked that the finite-size effects and the finite- 
temperature effects are much smaller than the symbol size in all of our results. Note how the 
QMC data follows the scaled curve much better than the strong-coupling curve, although 
there are definitely differences between the two. The deviations between the QMC data 
and the scaling result are real and larger than the finite-size or finite-temperature effects. 
This simply reflects the fact that the scaling result is not an exact interpolation formula 
for the momentum distribution. As expected, the momentum distribution is peaked at zero 
momentum, and as one approaches the critical point at x = 0.10224, the peak becomes 
sharper. One can also see that the truncated third-order expansion is not too accurate. As 
we already mentioned above, the curvature for momenta near the zone boundary has the 
wrong sign even for quite small hopping. It also underestimates the size of the peak at 
zero momentum, and this gets worse as we approach the critical point. Nevertheless, the 
strong-coupling expansion is quite accurate for small enough hopping, and the fact that it 
agrees essentially exactly with both the scaled results and the QMC simulations, provides 
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FIG. 8: (Color online) Momentum distribution function for the three-dimensional case with x = 
0.09 as a function of the band energy ek- Once again, the QMC data agrees better with the scaling 
theory than it does with the strong-coupling results or the scaled RPA, although deviations can 
still be seen in the data. Note that the scaled RPA works better than the truncated strong-coupling 
expansion. 

an independent check that all of these different approaches are working to high precision. 

We now try the phenomenological approach on the three-dimensional phase diagram. 
Here we have some uncertainty in how to proceed, because the scaling polynomial has 
freedom in our ability to vary the c' coefficient. We can either modify the scaling polynomial 
to represent the changes in c', or we can fix c' at a specific value and proceed from there. 
It turns out that we get better results if we fix c' = when calculating the phase diagram 
(especially for the two-dimensional case below). So we adopt that as our procedure (note 
we do not also set e' = 0, because that would produce the wrong critical point for this 
phenomenological approach). The result for the Mott phase lobes is 



1 1 2 ^ i - - 8.81514x3 ^ ^ 

= --x--x^ + x^± , ^ ^ . 100 

±2 2 VI - 12x + 20x2 + 16.67387x3 



U 

These results are plotted versus the QMC calculations [ITH in Fig. [9l One can see that while 

n Q 

the truncated strong coupling expansion [8|, |9| does not agree so well with the QMC data 
near the critical point, the agreement of the scaled curves is excellent. 



33 



ZD 
\ 



0.8 



0.6 



0.4 



0.2 - 



"1 ' 1 ' 1 ' r 

— Strong coupling 

Scaled 

• QMC 




0.02 0.04 0.06 0.08 
x=3t/U 



FIG. 9: (Color online) Phase diagram of the three-dimensional Bose Hubbard model. Note how 
the truncated strong-coupling expansion does not agree so well with the QMC data but the 
scaled results nearly fit the Mott phase lobe perfectly. 

Next, we move on to two dimensions. Recall that 7s = 0.64 in this case. Going tlirough 
the same procedure outlined above produces the following solution for the coefficients in the 
scaling polynomial: a = 8/7, = 4.6875; b = -17/2-/^ + 9/27^ = -2.29492; c = 48/7, + 
2c737s = 75.0 + 1.04167c'; d = 33/7, - 51/27^ + 9/27^ = 6.47278; and e = -256/7, + 
144/7^ - 2cV97, + 2cV7f + 2eV37, = -48.4375 + 4.53559c' + 1.04167e'. Once again, c' and 
e' are as yet undetermined. We find e' by requiring the critical point at k = to occur at 
the QMC and strong-coupling critical point Xc = 0.11948 
when 



16|. The critical point is found 



9.375xc + (9.57032 + 0.260418c')x^ 
(-27.56349 - 2.26780c' - 0.260418e': 



(101) 



X" 



0. 



(If we set c' = e' = 0, then the critical point would lie at Xc = 0.11579 which is a 3.2% error.) 
Substituting in Xc = 0.11948, then yields e' = —68.7054 — 0.338706c'. The parameter c' is 
then determined by requiring the integral of over all momentum to equal one. We find 
that c' ranges from approximately —115 at x = to c' ~ —224 at x = 0.119, but for values 
of X larger than about 0.1169, there is no value of c' that gives the total particle density 
to be exactly one — the error is about 1.5% at a; = 0.119 when we choose the best fit c'. A 
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FIG. 10: (Color online) Momentum distribution function in two dimensions with x = 0.05. We 
plot the strong-coupling expansion against the scaling theory results, the scaled RPA, and QMC 
simulations. Note how the QMC results agree much better with the scaled results and do not show 
the change in curvature near = 1- In addition, the scaled RPA doesn't work as well here as it 
did in three dimensions. 

simple fit of c'{x) is 



We compare our analytic expressions to QMC data in two dimensions on a 48 x 48 lattice 
with T = 0.05. In Fig. [10], we plot a case far from the critical point with x = 0.05. The scaling 
curve and the truncated strong-coupling expansion are both quite close to each other here, 
but one can see how the curvature has changed in the strong-coupling expansion but not in 
the data nor in the scaled curve. One also can see systematically that the QMC data agrees 
better with the scaled curve than the strong-coupling expansion. Moving on to a point much 
closer to the critical point at x = 0.1, we show the same plots in Fig. [TT]with the QMC data 
on a 48 X 48 lattice with T = 0.00625. Here, one can see a much more dramatic difference 
between the truncated strong-coupling results and the scaled results. While there definitely 
are some minor discrepancies with the QMC data and the scaled results, the agreement is, 
in general, outstanding. Note that we plot the momentum distribution versus instead of 
k, because in the strong coupling expansion all momentum dependence is summarized in ek 




(102) 
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FIG. 11: (Color online) Momentum distribution function in two dimensions with x = 0.1, which 
is close to the critical point. We plot the strong-coupling expansion against the scaling theory 
results, the scaled RPA, and the QMC simulations. Note how the QMC results agree much better 
with the scaling theory results. 

through third order, so there is limited other momentum dependence. For the QMC data, 
we average the small number of degenerate energy values. 

We finally try the phenomenological fit to the phase diagram by using the scaling poly- 
nomial in the power law and forcing the third-order strong coupling expansion to agree with 
the phenomenological scaling ansatz. Once again, we set c' = when we do this, because 
the agreement is significantly worse with different c' values. Because c' assumes much larger 
values in two dimensions in order to get the right integrated weight in the momentum dis- 
tribution, this is a significant assumption we are making, but as seen in the final results, the 
assumption seems reasonable because the agreement is quite good. 

Following an identical procedure to what was done in the three-dimensional case (with 
c' set equal to zero), we find 



U 



± 



1 
2 



X 



3 o 3 o 



(103) 



± 



\ + 0.14063X - 0.21460x2 - 3.87043x3 



[1 - 9.375X + 9.5704x2 + 9.6757x3]0-67 ' 
These results are plotted versus the QMC calculations |3| in Fig. [T2j Once again note that 
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FIG. 12: (Color online) Phase diagram of the two-dimensional Bose Hubbard model. Note how 
the truncated strong-coupling expansion [a, |9|] does not agree so well with the QMC data but 
the scaled results nearly fit the Mott phase lobe perfectly. 

while the truncated strong-couphng expansion does not agree so well with the QMC data 
near the critical point, the scaled curves lie essentially on top of the QMC data. 

The one-dimensional case is different from higher dimensions because the scaling behavior 
is not power law, but instead is the Kosterlitz-Thouless form of the two dimensional XY 
model. Hence, we modify our scaling ansatz to 



1 

X exp 



n 



-W 



w 



(104) 



which replaces the power law divergence by the appropriate exponential divergence. Because 
the exponent r] = 0.25 for the two-dimensional XY model, we have that W = 0.75W, with 
W the parameter in the Kosterlitz-Thouless fit to the one- dimensional Mott phase diagram. 

n 

Using the data of Elstner and Monien [Ij], we fit the gap function A(x) to the Kosterlitz- 
Thouless form 

2 _ A + Bx + Cx^ + Dx^ + Ex^ + Fx^ + Gx^ 
^^^^ ~ 1 + Hx + Ix^ + Jx3 + Kx^ + Lx^ + Mx^ + Nx'^ ' 

by using a Pade approximant for the pole that develops in the square of the logarithm 
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of the gap function. Note that one needs to do the Pade approximant for the square of 
the logarithm of the power series in order to obtain a robust fit [instead of doing a series 
or Pade approximation for A(x) first and then taking the square of the logarithm of the 
resulting series or Pade approximant]. The critical point is Xc = 0.29981 and the parameter 
W becomes W = 1.7241 or W = 1.2931. 
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FIG. 13: (Color online) Momentum distribution function in one dimension with x = 0.1, which 
is far from the critical point. We plot the strong-coupling expansion against the scaling theory 
results, the scaled RPA, and the DMRG calculations. Note how the DMRG results agree much 
better with the scaling theory results than the truncated expansion or the scaled RPA. 
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Now we solve for the coefficients in the scaling form just as we did in higher dimensions. 
First we ensure that the power-series expansion of the scaling form reproduces the strong- 
coupling expansion through the third order in x, then we ensure that the denominator of 
the square root in the exponential diverges at Xc- These two conditions yield a = 4.6400, 
b = 3.0006, c = 37.1201 + 1.0311c', d = 9.4879, e = 64.2632 + 6.8329c', and e' = -9.4630 - 
0.3190c'. The coefficient c' is adjusted to guarantee that the integral of the momentum 
distribution over all momentum is equal to one. We find that c' ~ —7.92 — 15.16x in order 
to satisfy the sum rule. 

We compare the scaled strong-coupling perturbation theory to the numerical calculations 
in one dimension from the density matrix renormalization group (DMRG) approach 37] 
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FIG. 14: (Color online) Momentum distribution function in one dimension with x = 0.2, which 
is two-thirds of the way to the critical point. We plot the strong-coupling expansion against the 
scaling theory results, the scaled RPA, and the DMRG calculations. Note how the DMRG results 
agree much better with the scaling theory results than the truncated expansion or the scaled RPA, 
but one can see that the scaling approach is beginning to fail. 

(provided to us by C. KoUath). Those calculations are essentially exact except for finite-size 
effects which become more important as we approach the critical point at x = 0.29981. 
In Fig. [ini we compare the different approximations to the DMRG calculations. One can 
immediately see that although the truncated expansion has a nonmonotonic dependence on 
ek, the scaled approach essentially agrees exactly with the DMRG calculations. 

Next, we compare the different approximate results to the DMRG calculations for x = 0.2 
in Fig. [HI Here we see that while the scaled results still agree well with the DMRG results 
near k = 0, the agreement is not so good throughout the Brillouin zone, and it is clear that 
the approximation is becoming inadequate. When we compare results for large values of 
the hopping, such as x = 0.25, the scaled results become negative over about half of the 
Brillouin zone, which is unphysical. 

We do not go through the phenomenological exercise of comparing our results to the 
phase diagram in the one-dimensional case as we did previously for higher dimensions. This 
is primarily because we could see the approximate scaled results were breaking down around 



39 



X ~ 0.2, so it is unlikely that a phenomenological approach for the phase diagram would be 
accurate in this case. In general, the strong coupling approach is more accurate in higher 
rather than lower dimensions. 

IV. CONCLUSIONS 

In this work, we have shown how one can generalize strong coupling perturbation theory 
from an expansion for the many-body energy levels, or for different ground-state correlation 
functions, to a direct expansion for the many-body Green's function at finite temperature. 
Here, we focused on applying the expansion to the problem of determining the momentum 
distribution in the bulk for the Bose Hubbard model within the Mott-insulating phase. By 
applying a scaling ansatz, that was motivated by recent work on the RPA, we are able to 
find accurate analytic expressions for the momentum distributions that hold nearly up to 
the critical point in two and three dimensions (the results for one-dimension are not quite as 
good). In addition, we showed how one can apply the results for the momentum distribution 
function to create a phenomenological theory for the Mott phase lobes. Comparing these 
results to QMC simulations showed excellent agreement in two and three dimensions. 

The strong coupling formalism as developed here can be used, as we have indicated, 
to obtain a strong-coupling expansion for the self-energy, and to include inhomogeneous 
features like a harmonic trap or disorder potential, and the effects of thermal excitations. It 
can also be readily adapted to nonequilibrium cases such as moving the origin of the trap 
or modulating the optical lattice depth for Bragg spectroscopy. The quantum Monte Carlo 
approach can be generalized to calculate dispersion relations, densities of states, and real 
time dynamics. We intend to examine those problems in the future. 

Acknowledgment s 

J. K. F. and H. R. K. acknowledge support under ARO Grant W911NF0710576 with 
funds from the DARPA OLE Program. H.R.K also acknowledges support from DST (India) 
as a J. C. Bose Fellow. Part of this work was completed during a stay at the Aspen Center for 
Physics. N. T. acknowledges support under ARO grant number W911NF0810338 with funds 
from the DARPA OLE Program. The quantum Monte carlo simulations were carried out 



40 



at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo. Y. K. 
and N. K. acknowledge support under MEXT Grant-in- Aid for Scientific Research (B) (No. 
19340109) and a Grant-in-Aid for JSPS Fellows We also acknowledge useful discussions with 
D. Arovas, A. Auerbach, C. Kollath, H. Monien, W. Phillips, J. Porto, R. R. P. Singh and 
I. Spielmann. DMRG data for the momentum distribution in one dimension were provided 
by C. Kollath. The QMC data for the phase diagram in two and three dimensions were 
provided by B. Capogrosso-Sansone. 

APPENDIX A: IMAGINARY TIME INTEGRALS NEEDED FOR THE 
STRONG-COUPLING EXPANSION 

Consider the i-time ordered integral 

o Jo Jo Jo 

It is easy to see that the sequence of functions satisfy the recursion relation: 

T (t-£ £ ■■■£)=[ r/r'p"(^~^'^^°o7 ,(t'-£ £ ■■■ £ ) (A2) 

Jo 

Taking the Laplace transform of both sides, it is straightforward to see that 

poo 

^[■^m{Ti £aoi ^Omi ' ' ' i^ai)iS] = / d-TC Im{Ti £aoi £ami ' ' ' 

^0 

POO PT 

Jo Jo 

POO POO 

Jo -5 + ^ao 

^[^m-liT] £ami ^Om-i ' ' ' i£ai)]s]. (A3) 



■5 ~l~ ^'qq 



Iterating this, and noting that /o(r;£^Q,J = e "^^"i/S, which implies that C[Io{r; £aj; s] = 
[3(s + £aj]~^, we find 

£=0,171 ^ 
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Taking the inverse Laplace transform yields 




'^+*°° ds e 



,TS 



n 



1 



(A5) 



2m 3 



e=o,m 



with 7 > max{Sao,Sa„^, ■ ■■ ,Sai), so that all the singularities of the integrand lie to the 
left of the integration contour in the complex s-plane. The integral is straightforwardly 
evaluated using the contour integration techniques. When all the energies Sao,Sam, ■ ■ ■ , £ai 
are distinct, we get one contribution from each of the m + 1 simple poles of the integrand in 
Eq. flASp . leading to Eq. fl5T]) . If one and only one pair of energies are equal, say, £a^ = Sap, 
then the integrand of Eq. (lASp has m — 1 simple poles and one double pole, and we get 
Eq. ( l52l) . One can similarly extend the results to other cases, corresponding to two double 
poles, or one triple pole, etc. 
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APPENDIX B: FINAL RESULTS FOR THE THIRD-ORDER EXPANSION 
TERMS 



Explicit forms for the third-order coefficients in the strong-couphng expansion are pre- 
sented here (for brevity only in the second form, as discussed in Sec. II): 
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Similarly, we get 



^jj2jj' ~ ^jh^j2Ajj' Pj,n Pj2,n2 Pj',n' 



+ 
+ 



71,712 ,n' 



n{n2 + l)n{n' + 1) 1 



+ 



(n + l)(n2 + l)(n + l)n' [n + l)n2(n + l){n' + 1) 



~'~ ^j',n')(^j2,n2 ^j' ,n')i^j,n + S''.'"''' ^S," ~'~ ^32,n2)^^j,n + ^32,^12)^^32,^12 ^j',n') 

, nn2{n + l){n' + 1) _ 1 1 ^ 



+ 
+ 



^^j,n ~^ ^j',7i') (*^j2i"2 ~^ ^3,n) ~^ ^j',n') ^^32,''^2 ~^ ^j,''^ 

{n + l)?7,2(n + 2)?T,' 

(^j'^Ti + S'2,n,2 + + S-',7i')(^j'^n + ^ j ' ,n') i^t,n + S''."'-' 

n(n2 + l)(n - l){n' + 1) ^ 

(^j,n ^32,^2 ~^ ^j,n—l ~^ ^j' ,n')i.^j,n ~^ ^j2,n2)^^j2,n2 ~^ 



. nn2(n + l)(n' + 1) n(n2 + + l)n' 



(^jr,71 ~^ ^j',7l') (^J2,712 ~^ ^J,71.) (^J2,712 ~^ ^j',n') ^^i,71 ~^ ^j2,12 ) (^i2 ,"-2 ~^ ^i',7l' ) (^J,71 ~^ ^j',7l') 

n + l)n2(n + 2)n' 



+ 
+ 
+ 



(^J,7l ~^ ^J2,712 ~^ ^iiJ^+l ~^ ^i',7l' ) (^i,71 ~^ ^^'2,712) (^i, 71 ~^ ^j',n' ) 

n(n2 + l)(n-l)(n^ + l) ^ 

^^j,n ~^ ^32,1^2 ~^ ^i,"— 1 ~^ ^ j' ,n')^^j 2,n-2 ~^ ~^ ^j',n') 

^ n(n2 + l)(n- + 1) 

(n + l)n2(n + 2)n' 

^^j,n ~^ ^32,^2 ~^ ^iiJ^+i ~^ ^i',7i') (^i2i"2 ~^ 

{n + 1)(?T.2 + (n + l)n2n(n' + 1) 



(^J2,n2 ~^ ^ j ,n) i.^ j2 ,n2 ~^ ^j',n')^^j,n ~^ ^j',n') ~^ ^ 32,^12^^^ 32,^12 ~^ ^j',n')^^j,n ~^ ^j',n') 

^ r2(n2 + l)(n- + 1) 

(^■,71 + ^j2,n2 + ^/.nO'-S'.n ~'~ ^i',71' ) (^'.n + ^j',n') 

(n + l)n2(n + 2)n' 



(^J,7l ~^ ^i2,7l2 ^ ^ ^j',7l' ) (^J,71 ~^ ^^2,712) (^i2,l!.2 ^ ^Ji"^ 

(n + l)(n2 + l)W 1 1 

(^J2,n2 ~^ ^J,7l)(^jr',71 ~^ ^j',n') ^^^21^2 ~^ ^ii"^ ^^ii" ~^ ^j',n') 



nn2n{n' + 1) ?t-(?t-2 + 1)^^' 



(^7,71 ^7', n') (^12,712 ~^ ^I'.n') (^7,71 ^7' 71') (^7,71 ~^ ^72,712) (^i2>"2 ~^ ^3,1^)^^32,^2 ~^ ^i'.n') 



Ji" ' 3 ,1^ 32,n2 ' ,n' J \'-j,n ' '-j',n'J \ 3,""- ' 32,n2' \"32,n2 ' 3,'"-l \"32,n2 ' 

^ (n + l)n2(n + l)n^ 1 ^— )] } 

(^j,7i ~^ ^i2>™2^ ^^ij*^ ~^ ^^3,1^ ~^ ^32,1^2} ^^ii" ~^ ^i'l"'^ 
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X 



^ n{n' + + 1) 

(n + l)n'(n + l)n' 



\ j,n ' j',n'. 



nn'{n + + 1) 



(/5 



{n+l){n' - l){n + 2)n' 
r2(ra' + 2)(n- l)(r2' + 1) ^ 



nn'{n + l)(n' + 1) 



+ l)(n + 2)n' 

n(n' + l)(n-l)(n' + 2) ^ 

(^jr> + ^i',n' + ^j,n-l + ^j'',n'+l)(^j'',n' "I" ^j,n)i^j,n + ^j'',n') 

^ n(n' + 2)(n- l)(n' + 1) 

(ri + l)(n + 2)n' 



{n + l)r;,'?7,(r;,' + 1) 



^ n(n' + 2)(n-l)(n' + l) 

(n + l)(n'-l)(n + 2)n' 



{n + l)(n' + l)nr;,' 
n{n' + l)n{n' + 1) 



(n + l)n'(n + l)n' 



iP 



-)] }• 
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The symmetry of the various terms in, and the term by term correspondence between, 
Eqs. flBH - |B4|) above are noteworthy. The above results are sufficient for the purposes of 
this paper, where we discuss only bipartite lattices (specifically, hypercubic lattices in d 
dimensions with nearest-neighbor hopping only). 
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